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We have numerically and analytically studied ac+dc driven 
Josephson-junction arrays with a single vortex or with a sin- 
gle vortex- antivortex pair present. We find single- vortex steps 
in the voltage versus current characteristics (TV) of the ar- 
ray. They correspond microscopically to a single vortex phase- 
locked to move a fixed number of plaquettes per period of the 
ac driving current. In underdamped arrays we find vortex 
motion period doubling on the steps. We observe subhar- 
monic steps in both underdamped and overdamped arrays. 
We successfully compare these results with a phenomenolog- 
ical model of vortex motion with a nonlinear viscosity. The 
I-V of an array with a vortex-antivortex pair displays frac- 
tional voltage steps. A possible connection of these results to 
present day experiments is also discussed. 

PACS numbers: 74.50.+r, 74.60.Ge, 74.60.Jg 



I. INTRODUCTION 

The presence of giant Shapiro steps and giant frac- 
tional Shapiro steps in the I-V characteristics of 2-D 
Josephson-junction arrays has attracted significant at- 
tention recently 0-0] . These 2-D arrays may be of use as 
a source of coherent microwave radiation g. In a sepa- 
rate context the flux-flow dynamics of vortices has been 
studied P-pT|. The reported experimental observation 
of ballistic vortex motion has also stimulated fur- 
ther theoretical and numerical investigations |l3|-[l7[| of 
the mass and friction of a vortex in an array. Until now 
these numerical investigations have focused on dc driven 
vortices. In this work we perform numerical simulations 
on Josephson-junction arrays, with only one vortex or 
with a vortex-antivortex pair present in it, driven by a 
time-dependent current i{t) = idc iacCos{27rut). We 
calculate the voltage V versus idc characteristics {I-V). 
We find harmonic and subharmonic single- vortex voltage 
steps and analyze the underlying phase- locked vortex mo- 
tion. A vortex-antivortex pair separated by a distance 
Ax along the direction of the injected external current 
phase-locks on to fractional voltage steps. 
The arrays are 2-D lattices of superconducting islands 
(sites) connected by Josephson junctions (bonds). The 
unit cells (plaquettes) of these lattices can be, for exam- 
ple, square or triangular. The vortices are represented by 



eddy-current patterns about a plaquette. Here we con- 



sider the classical regime defined by Ej Ec 



72C, 



where Ej is the Josephson coupling energy and Ec the 
charging energy of two islands, e the electron charge, and 
C the capacitance of a junction. In this regime quantum 
fluctuations are neglected, leaving the phases 6{r) of the 
Ginzburg-Landau order parameter on the islands as the 
only dynamical variables. 

In this case the array is well-modeled by the Resistively 
Capacitively Shunted Junction (RCSJ) model, defined by 
the total bond current i{r^r') between nearest neighbor 
sites r and r' 



i{ry)=M{ry) 

+ 0{r,r') +sin[6>(r,rO 



2nA{ry)l (1) 



plus Kirchhoff 's current conservation conditions at each 
site. Here the dots represent time derivatives. The three 
contributions to i(r^r') are the displacement, the dis- 
sipative and the superconducting currents, respectively. 
The phase difference across a junction is ^(r, r^) = 0{r) — 
Oy). The currents are expressed in units of the junction 
critical current Ic] time is measured in units of the char- 
acteristic time l/uoc = h/{2eRnIc)^ and Pc = (^c/^p)^ is 
the Stewart-McCumber parameter [Q, with the plasma 
frequency cOp defined as cOp = 2eIc/hC^ and Rn is the 
junction's normal state resistance. The bond frustration 
variable A(r, r^) is defined as the line integral of the vec- 
tor potential A: 



<P0 Jr 



A-dl, 



(2) 



with the elementary quantum of flux = h/2e. The 
frustration parameter / measures the average flux pierc- 
ing a plaquette, measured in units of (j)o. 
The motivation for this paper is twofold: to study the dy- 
namics of a few vortices in an array, and to see whether 
the results can be generalized in order to explain dynam- 
ical non-equilibrium states, like the axisymmetric coher- 
ent vortex state Here we deal mainly with the anal- 
ysis of single- vortex voltage steps. The motion of a vor- 
tex produces a Faraday voltage across the array. In this 
paper we find three types of new steps. First we find 
single- vortex voltage steps. The voltage V on these steps 



1 



is proportional to an integer multiple n of the frequency 
of the ac drive: 



V 



nhu 
~2e 



n = 0, 1, 2,-- 



(3) 



In the following we will also consider natural units h/2e = 
1. Then the voltage is normalized such that it corre- 
sponds to 27r times the number of jumps between plaque- 
ttes of the vortex per time unit. This differs from another 
often used normalization by a factor Ny^ the number of 
junctions perpendicular to the direction of the current 
injection. On such a step the motion of the single vortex 
is phase-locked to move an integer number of plaquettes 
per period ^. Next we find subharmonic single- vortex 
steps in the /- Vs of both overdamped and underdamped 
arrays. On these steps the voltage is: 



^ n hv 
m 2e 



n= 0, 1, 2, 
m= 1, 2, • • 



(4) 



The dynamics on these steps corresponds to the vortex 
moving n plaquettes in m periods. 

Finally we simulate a vortex- ant ivortex pair and obtain 
the following steps: 



V 



nNy hv 
m 2e 



n= 0, 1, 2,- 
m= 1, 2, •• 



(5) 



Each vortex moves nNy/ 2 plaquettes every m periods. 
These new steps should be contrasted with other types 
of steps that arise in Josephson-j unction arrays. In single 



junctions, harmonic (V 



' 2e 



) voltage steps are present 



in both experimentally and numerically obtained /-F's 
pT| . These steps are called single-junction Shapiro steps. 
On these steps the single junction has n phase slips per 
period ^. In simulations of and experiments on ac+dc 
driven arrays so-called giant Shapiro steps have been ob- 
served On these steps the NxNy individual junc- 
tions along the direction of the external current are all 
phase-locked on the same single junction Shapiro step, 
and one obtains 



V 



N,Ny 



nhv 
~2^ 



0, 1, 2,- 



The harmonic and subharmonic single-vortex voltage 
steps are therefore small in comparison to the giant 
Shapiro steps. Subharmonic giant Shapiro steps are ob- 
served experimentally and in simulations [|U^] . A par- 
ticular example of these half-integer steps is found in the 
axisymmetric coherent vortex states (ACVS). These sta- 
tionary states correspond to an oscillating pattern of vor- 
tex and ant ivortex streets, which arrange themselves at 
a well-defined angle with respect to the current direction 



In studies of and experiments on ac+dc driven arrays, 
with an average rational flux piercing through each 
plaquette, giant fractional Shapiro steps were observed 
On these steps the voltage is: 



NxNy nhu 

y — 

a 2e 



0, 1, 2, 



An explanation for such steps has been proposed and ver- 
ified in simulations : an / = | array has q degenerate 
ground states, consisting of vortex lattices, carried into 
each other by translations. In one period of the driv- 
ing current the vortex lattice moves from one degenerate 
state to the next. After q periods every vortex has moved 
through the whole array, generating the observed voltage. 
This shows that there are two kinds of Shapiro steps in 
Josephson-j unction arrays. Those based on the coherent 
phase slips of all the individual junctions in the array 
(Giant Shapiro steps) and those involving coherent oscil- 
latory vortex motions (giant fractional Shapiro steps and 
ACVS). 
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FIG. 1. Array geometry used in the simulations, illustrated 
with an 8 X 8 array. Junctions are denoted as crossed bonds. 
Periodic boundary conditions are imposed in the ^/-direction, 
while the current bias ib is applied along the x-direction. 

The approach of our study is to have only one vortex, 
which makes it is possible to separate the effect generic 
of vortex motion from the effects of interaction between 
them. One can systematically study the underlying mi- 
croscopic dynamics of the phase-locked vortex motion. 
Then by considering a vortex- ant ivortex pair one can 
study the effect of interaction in its most simple form on 
the phase-locked states. The vortex dynamics on single- 
vortex voltage steps has a number of new and interest- 
ing features, e.g. we observe vortex motion with period 
doubling in underdamped arrays. The underlying micro- 
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scopic vortex motion repeats itself only after 2, 4, 8, or 
even 16 periods ^, although the vortex still moves, on 
average, a fixed number of plaquettes per period. In this 
paper we discuss how these steps can observed experi- 
mentally in arrays with a low density of vortices. 
Previous authors have studied simple models of 

vortex motion. In these models the vortex is described 
as a point particle of mass M(/3c) experiencing a certain 
friction, and moving in a sinusoidal potential. In Ref. 
]T7[ we found that the I-V of a dc driven array in the 
vortex regime is well described by such a vortex equation 
of motion in terms of the vortex-coordinate y containing, 
instead of the usual linear viscous force, a nonlinear one 

M{Pe)y + ^ fsii^-^iyi + ^0 sin 2Try = % (6) 

In this equation A(/3c) and B{/3c) are phenomenologi- 
cal parameters. This analysis is similar in spirit to the 
description of long Josephson junctions in terms of the 
fiuxon coordinates An interesting question previ- 

ously left unaddressed, and considered in this paper, is 
whether a phenomenological model for the vortex coor- 
dinate can still reproduce the I-V of an ac+dc driven 
array. In this work we resolve this question and we com- 
pare the I-V of an array with the result of Eq. (^) with a 
time-dependent current i = idc 'i^ac cos 27ri/t, and using 
the parameters M{/3c)^ ^{Pc) and B{/3c) obtained for dc 
driven arrays. The results of Eq. (|^) and the simulations 
are in reasonable agreement over a broad range of values 
of the frequency and the amplitude of the ac drive. In 
other words the vortex experiences a nonlinear friction in 
an ac+dc driven array. The outline of the paper is as fol- 
lows. In section II we discuss the calculational algorithm 
to compute the I-V. In section III we present the I-V^s 
containing the harmonic and subharmonic single-vortex 
steps and discuss the possibility of experimental verifica- 
tion. We then study subharmonic single-vortex voltage 
steps and period doubled vortex motion using the micro- 
scopic current distribution in the array as a function of 
time. Finally we compare the result of Eq. (^) to those 
of the /-V's obtained from the Josephson-junction array 
simulations. In section IV we investigate the effects of in- 
teraction and discuss the vortex- ant ivortex voltage steps 
in the I-V. In section V we present our conclusions. 



II. CALCULATIONAL APPROACH 

We numerically solve the equations of motion for a two- 
dimensional array. We show the square lattice of Lx x Ly 
sites in Fig. |l|. The sites are connected through Joseph- 
son junctions, denoted by crosses. We use the RCSJ- 
model of Eq. ^ to relate the current i(r^r') through 
the Josephson junction to the phase difference 0{r^r'). 



We use periodic boundary conditions (PBC) in the y- 
direction, while the current is fed in and taken out along 
the x-direction. This set of coupled nonlinear differen- 
tial equations can be integrated efficiently using a fast 
Fourier transform algorithm |^,|2^ . The array consists of 
Nx X Ny plaquettes {Nx = Lx — I and Ny = Ly). The 
vorticity n{R) is defined as: 

27^n(i^) = 27r/ + ^ ((9(r , r') - 2iiA{r, r')) . (7) 
ViR) 

Here V{R) denotes an anti-clockwise sum around the 
plaquette R and the gauge- invariant phase difference 
0{r^r') — 27rA(r,r^) is taken between — tt and +7r. 
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FIG. 2. /-l/'s obtained from simulation of a 16 x 32 array 
with parameters, iac = 0.10, v — ^ and fic — ^ (continuous 
line); v — ^ and /3c = 5 (dotted line); and ^ — ^ and 
/3c = 20 (dashed-dotted line). The voltage in natural units is 
normalized by a; = 27rz^. The inset contains a close-up of the 
/3c = 20 curve, with an added I-V branch showing hysteresis. 

We are interested in the behavior of vortices in 
Josephson-junction arrays. Stable vortices can be ex- 
plicitly introduced in the initial phase-configuration by 
the application of a small frustration . The plaquette 
coordinate II with unit vorticity will be called the topo- 
logical vortex coordinate. The voltage V(i) is obtained 
from: 

E j^^{Lx-\.y)-my)\ (8) 

The time average of V(t) is related to the average vortex 
velocity ^ = The microscopic dynamics of the vor- 
tex motion is refiected in the eddy-current distribution 
C(R,i), 

C(il,t)= ^ i(r,r',t). (9) 

V(R) 
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Where V{R) is the anti-clockwise sum over bonds about 
a dual lattice site R. The vortex shows up as a local 
extremum in the eddy-current distribution. 
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FIG. 3. Period doubling on the n = 1 step. The results are 
obtained from a simulation of an 8 x 8 array, with parameters 
f3c = 25, = 1/50, iac = 0.10, and idc = 0.13. In (a) we 
plot the vortex position versus time. In (b) we plot the eddy 
currents versus time at three adjacent plaquettes of the middle 
column. The continuous, dotted and dot-dashed lines are 
respectively the first, second and third plaquette in the middle 
column of the array. The minima in these graphs indicate the 
position of the vortex. 



III. SINGLE- VORTEX-INDUCED VOLTAGE 
STEPS 

In this section we present the results of our simulations 
of a single vortex in a Josephson-j unction array. In the 
first subsection we discuss results for the I-V character- 
istics. In the next subsection a detailed description of 
the single-vortex voltage steps is given. We discuss the 
microscopic vortex motion, including its period doubling 
and the subharmonic single- vortex voltage steps. This is 
followed by a comparison of a simple model for vortex 
motion to the simulation results, and finally the discus- 
sion of possible experimental verification. 



plaquettes. 

The second branch of the I-V is generated by decreasing 
the bias current. As the initial phase configuration one 
uses the final phase configuration obtained in the upward 
sweep. At the retrapping current irt the average voltage 
returns to zero. In under damped arrays the retrapping 
current irt can be different from i^p. This hysteretic be- 
havior has often been interpreted as evidence for inertia 
of the vortex in the dc driven case. An example of hys- 
teresis is shown in the inset of Fig. ||. Focusing the dis- 
cussion on the upward current sweep we encounter the 
first plateau at V = 27tu. This is the first single-vortex 
step. The vortex moves on average one plaquette per 
period. There are two more steps visible at multiples 
of 27ru in Fig. ||. Above idc ^ 0.25 no steps are visible 
anymore. The step width has become smaller than the 
current grid size Aide = 0.01. The /3c = I-V exhibits 
a pronounced upward curvature. Hence the vortex vis- 
cosity is nonlinear [|l^. Between the n = 1 and n = 2 
step a small subharmonic step, at V/27rh' = |, is visible 
in Fig. ||. On this step the vortex moves three plaquettes 
every two periods. Using a smaller Aide = 0.002 one can 
even observe the | and | steps. 

Increasing the damping parameter Pc for given shifts 
the steps to higher values of the current bias idc- This 
is due to an increase in the vortex viscosity with Pc 
| ]T7[ . When (3c is changed from 5 to 20, the n = 3 step 
width is gradually reduced to below the current grid size 
Aide = 0.01. For u > 1.0 no integer steps are present 
any more. For these frequencies the velocity for which 
the vortex would phase lock on the first integer step lies 
above the maximum vortex velocity in the array. 
The width of a particular step varies in an oscillatory 
fashion as a function of j^. This is qualitatively similar 
to the step width behavior of a single junction, which 
varies as a Bessel function of iae- By adjusting iac one 
can make more steps visible. 



B. Single- vortex-induced voltage steps 



A. The I-V's 

The /-U's are generated by gradually increasing the 
bias current ide from starting with an initial configu- 
ration containing a vortex. We express the voltage in 
natural units ^ = 1 and normalized by the frequency 
27TU. The harmonic steps then occur at integer voltage 
values. 

For currents below the depinning current idp the voltage 
is zero (see Fig. |^), when it is averaged over enough pe- 
riods of the ac drive. The vortex deforms in response to 
the ac+dc drive, but stays in the same plaquette [p3| . 
Or, for low enough frequency u and large enough i^c, it 
can even oscillate back and forth over a finite number of 



We now turn to the microscopic vortex motion on the 
single- vortex voltage steps. We first consider the n = 
step, i.e. the response of a vortex that is pinned in one 
plaquette. One may expect the response to have the same 
period as the ac drive. In that case the quantity Vp^ the 
voltage averaged over one period of the drive, is constant 
and equal to zero. We find, however, that Vp can be non- 
zero on the n = step. This is repeated periodically in 
time. In table || we show the periodicity of Vp^ (in units 
of the driving period) as a function of ide for an 8 x 8 
array with parameters /3c = 25, v — 1/50, iae = 0.1. 
Next we focus on the n = 1 step in table |[ The vortex 
depins at ide ^ 0.12. For ide = 0.13 the voltage Vp alter- 
nates between two different values. Figs. || and || show 
the corresponding vortex motion. Fig. 0(a) shows the 
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topological vortex coordinate, defined in Eq. (|^), versus 
time. It shows that the time the vortex spends at one 
plaquette alternates between two values Ti and T2. The 
sum Ti + T2 is equal to In Fig. ||(b) we plot the 

eddy currents of three adjacent plaquettes versus time. 
In Fig. ^ we show snapshots of the spatial distribution of 
eddy currents. Frames and 10 are almost equivalent. A 
jump occurs between frames 2 and 3 and between frames 
8 and 9. The difference between Ti (frames 0-2) and T2 
(frames 3-8) is clearly visible. All the figures show that 
the actual vortex motion is periodic with twice the pe- 
riod of the driving force. 
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FIG. 4. Period doubling on the n = 1 step. Parameters are 
the same as in Fig. 3. We show snapshots of the eddy current 
distributions at different times. The time step between two 
snapshots is 10 (in units of l/uc). The eddy-current distribu- 
tions are smoothed out by interpolation. 

Subharmonic steps turn out to be a generic feature of 
the I-V^s of ac+dc driven arrays. They are observed for 
different values of iac, ^, Pc and system sizes. Many dif- 
ferent voltage values are possible. We studied a | step 
in detail. At idc = 0.11 we found an n = ^ step in a 
/^c = O7 ^ = 1/25, 8 X 32 array. On this step the generic 
behavior is as follows. Although the vortex deforms sig- 
nificantly during the first period, the topological vortex 
coordinate does not change. Only in the second period 
it jumps to the neighboring plaquette. We now describe 
a typical example. 

In Fig. H we show the eddy current versus time of the 
plaquette containing the vortex. Both curves have two 
minima: A and B. In minimum B the vortex has a dif- 
ferent shape (and hence feels a different potential) than 
in A. From this we can deduce the following scenario. 
During the first period the vortex deforms from config- 
uration A to 5, staying in the same plaquette. In the 



following period it jumps to the next plaquette into the 
adjacent minimum A. We also found a case {idc = 0.072, 
u = 1/50) in which the vortex briefly jumped to the next 
plaquette (and returned to the original plaquette) during 
the first period. 
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FIG. 5. Vortex motion on a ^ step for an 8 x 32 array, with 
parameters /3c = 0, = 1/25, iac = 0.10, and idc = 0.11. 
We plot the eddy current versus time, as in Fig. 3(b), for 
two adjacent plaquettes (denoted by the continuous and the 
dot-dashed line respectively). The minima A and B in this 
graph indicate the position of the vortex. 

C. Comparison to a phenomenological model of 
vortex motion 

In this subsection we compare the results of the numer- 
ical simulations to the results of a simple model for the 
vortex motion. The motion of a single vortex in an ar- 
ray can be modeled by that of a point particle with mass 
M that, driven by a (Lorentz) force z, moves through a 
sinusoidal pinning potential and experiences a linear vis- 
cous damping force with constant viscosity coefficient 77 
flH • The vortex mass M can be calculated by equating 
the electromagnetic energy stored in the array to a vor- 
tex kinetic energy ^My'^. The friction can be determined 
by equating the dissipated energy to r]y'^. This results in 
the following equation of motion for position y{t), gener- 
alized to include a time dependent driving force i. 

TTpci) -\- Try -\- io sin 27ry = i^c + iac cos 27rut. (10) 

For a square array ^ 0-10 Q. Every time the particle 
moves one plaquette (?/ ^ ^ + 1) an integrated 27r voltage 
pulse is generated across the array in the x-direction. The 
average of 27ry is then the dimensionless average voltage 
measured across an array. 

We recently introduced a modified vortex equation of mo- 
tion for dc driven vortices, including a nonlinear vis- 
cosity 

A(6 )y 

^if^c)y + , ^| ., + ^0 sin 27ry = idc + iac cos 27nyt. 

(11) 
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The parameter values for overdamped case are M{/3c) = 
0, A ^ 2.7, B ^ 1.8, io = 0.1 0. 

Equations (|lO|) and ( |ll| ) are similar to the equations de- 
scribing single Josephson junctions. To connect to the 
single- Josephson-j unction literature, replace y by a phase 
= 27Ty^ divide by io and absorb the coefficient in front 
of the ffist derivative of ^ in a new time unit: 



sin ^ = idc + iac cos l^r. 



Where /3 = M{Pc) 



io27r 



B 



(12) 



and 



— Idc 



To allow for a more quantitative comparison we use the 
procedure described in Ref. to find the parameters 
A and B in Eq. ([ll|). That is we fit the I-V curve in the 
dc driven case to the form 



dc 



Given these parameters one can numerically calculate the 
I-V predicted by Eq. (pT[). In Fig. || we compare these 
results for one specific iac^ ^ and /3c to the I-V obtained 
from simulations of the full array. The results of Eq. ( pj] ) 
are in good agreement with the simulation for larger cur- 
rents idc > 0-3. For this current regime the steps are too 
small with respect to the current grid Ai^c = 0.01 to be 
visible. 
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FIG. 6. Comparison of results the model Eq. (12) (dotted 
line, parameters: Q — 7r/5, f3 — 0.4, and iac — 1-0) and 
(full line, parameters: Q = 0.617, B = 0.0424, /3 = 0, and 
iac — 1.0) to the simulations of the full 16 x 32 array (circles, 
^,/3c = 2, and w = 0.10). 



For lower currents the agreement is less satisfactory. 
The steps in the simulations and in Eq. ( pj] ) do overlap 
to a large extent, but the lower edge of the steps is under- 
estimated by Eq. (pj]), especially at depinning. However 
one finds large deviations when comparing the simulation 



data to Eq. (p^oj), the model with linear viscosity. 
The phenomenological vortex mass for pc < 35 is 
M{pc) = 0, as found from the dc I-V curves in Ref. [Q. 
In the dc-hac case, however, we find small hysteresis loops 
in the I-V. One could interpret this as an indication for 
the presence of a non-zero M{/3c). The hysteretic behav- 
ior of Eq. (|l^) is complicated. Using different M 7^ 
with the computed values of A{Pc) and B{Pc) did not 
yield better agreement with the simulation results. 



D. Experimental verification 

In order to see how the single vortex phase-locking 
mechanism would manifest itself in experimentally ac- 
cessible conditions, we now reinstate physical dimensions 
using experimental parameters [^. Typical values for 
the parameters are Ic = 0.01 — 2.0 /iA, IcRn = 300 jaV, 
oJc = ^ 1 GHz and pc = 10-100. The single- 

vortex voltage steps would occur for i^c = 0.10 — 0.30 Ic 
and for frequencies u < 0.10-cc;c ~ 100 Mhz, and < 50. 



Experiments are conducted at a finite temperature. Volt- 
age steps are observed as a reduction in the differential 
resistance while at zero temperature ^ = on a 
step. This reduction should be large enough, and ex- 
tend over a sufficiently large current range, in order to 
be measurable. The temperature is expressed in units of 
To = 2^ ^ 2 X 10^ Ic. We have performed a T 7^ 



simulation for one particular case: a 16 x 16 array with 
parameters = /^c = 0, and iac = 0.10. The simu- 
lation was done using the algorithm introduced in Ref. 
| p5[ , as extended for arrays [^. We found that up to 
T = 0.004 To the n = 1 step was clearly visible. For 
T = 0.008 To a reduction in ^ was hardly discernible. 
This puts an upper bound on the appropriate tempera- 
tures, varying from T = 1 to 200 mK with the value of 

Ic- 

To study the single- vortex phase-locking phenomena ex- 
perimentally, one can use relatively small samples, to 
make sure that only a small absolute number of vortices 
is present at any time, or a large array with a low vortex 
density. In the former type of experiments the interaction 
of the vortex with the boundaries, or equivalent ly with 
image antivortices, has to be taken into account. When 
moving to a boundary, a vortex may either escape from 
the array or reffect as an antivortex jioj . In the parame- 
ter regime studied in this paper, vortices are entering the 
array and leaving it again at the opposite boundary. In 
practice the small number of vortices present will there- 
fore ffuctuate slightly. This may effect the way in which 
phase-locking is established as compared to our simula- 
tion in which only one vortex is present at any time. To 
estimate this effect we have simulated a finite array in 
a small magnetic field and also found steps in the I-V's 
due to phase-locking with the ac drive at voltages slightly 
different from those mentioned above. 
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We also have performed simulations of triangular ar- 
rays, often used experimentally, and find qualitatively the 
same scenario for the occurrence of single- vortex voltage 
steps, as reported in section IIIA. 



1.0 




FIG. 7. I-V of a vortex and an antivortex separated by a 
distance Ax = 4, 5, 7, 8, and 10 in an overdamped dc driven 
32 X 32 array. The annihilation current is the current at which 
the average voltage drops to zero in the I-V of the downward 
current sweep. 



IV. VORTEX- ANTIVORTEX PAIR-INDUCED 
STEPS 

In this section we consider the effect of vortex- 
antivortex interaction on the single- vortex voltage steps. 
In the first subsection we show the fractional voltage 
steps in the I-V due to the vortex- antivortex interaction, 
and obtain these steps from a phenomenological model 
for vortex motion including the logarithmic interaction. 
In the second subsection we briefly discuss the analogous 
behavior of excess and missing vortices in a checkerboard 
groundstate of an / = | array. 

A. Vortex-ant i vortex pair in an / = array 

We calculate the zero temperature I-V^s of an array 
containing a vortex- antivortex pair. There is no applied 
magnetic field, and the pair is included in the initial con- 
figuration by construction In order to show how 
one may obtain meaningful results starting from such 
a metastable configuration, let us first discuss the case 
of only a dc drive. When the applied current is zero, 
the only mechanism that may prevent the annihilation 
of the pair is the pinning of the lattice. Pinning will pre- 
vent annihilation if the mutual separation is at least 8 
lattice constants. This minimum distance is in approx- 
imate agreement with the value of 0.10 (in units of 



the junction critical current) for the maximum pinning 
force on a vortex in an infinite system. However, when 
the current is non-zero, and the vortices move in different 
rows (in opposite directions), annihilation may be absent 
even when the perpendicular distance between the rows 
is less than 8 lattice constants, depending on the magni- 
tude of the current. We have studied this stabilization 
by the driving force quantitatively by performing down- 
ward current sweeps starting from a number of initial 
configurations containing pairs with different separation 
distances Ax. The /-Vs are shown in Fig. |^. When the 
pair is unstable, and thus annihilated, the average volt- 
age obtained from the simulations sharply drops to zero. 
This annihilation current strongly depends on the sepa- 
ration distance of the pair in the initial configuration. A 
higher current is needed to stabilize a smaller pair. 
Next we consider the behavior of an ac+dc driven pair. 
Since two vortices move in opposite directions under the 
influence of the current, their mutual distance will change 
in time. Phase-locking can therefore not be established 
for a voltage corresponding to one jump per period, since 
in the next period the environment has changed, and 
thereby the interaction strength. We nevertheless find 
steps in the I-V of these systems. Once the vortex and 
the antivortex have traversed Ly/2 plaquettes, they cross 
each other again because we have periodic boundary con- 
ditions. Phase-locking is established if at that point the 
phases are the same again. If phase-lock is established 
over n traversals, and the motion is periodic with period 
m (in units of the external frequency u), then the ob- 
served voltage is: V = nLy/m27ru. We can find each of 
these steps in the I-V in Fig. H, although for higher n the 
step width decreases rapidly. The step width is typically 
smah (- 0.01). 

When the vortex and antivortex are separated by a large 
distance, one would expect to see single-vortex voltage 
steps again. This does not occur for the system size we 
considered due to the long range of the interaction. In 
section IIIB we considered models for a single vortex, 
here we consider a model for a pair. We model the vortex- 
ant ivortex interaction as being logarithmic. In addition 
we have to take into account the periodic boundary con- 
dition along the ^/-direction: the vortex can see an an- 
tivortex behind it and in front of it (and vice versa). 
We denote the position of the vortex by and that of 
the antivortex by y2. Let the coordinates take the val- 
ues < i/i < Ly. The images are then at i/i Ly. We 
only take the direct interaction and the interaction of the 
nearest image charge into account (i.e. at a distance less 
then Ly). The only stable current-driven pairs in the 
simulation are the ones that move at a fixed separation 
Ax, we therefore fixed the separation Ax in this model. 
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This leads to the fohowing equation of motion: 



^(^i)^i = - ^dsin(27r?/i)- 



yi - y2 



Ax' + (^1 - 



Ly - \yi - y2\ 



V{y2)y2 = -i{t) - idSm{27Ty2)- 



yi-y2 



- \yi - y2\ 

'Ax'^{Ly-\y,-y2\r 

Here e is the sign of yi — ^2- We have computed the I- 
Vs using these equations and we indeed found steps at 
V = nLy/m27ru. To capture the mechanism for the pair 
steps in the model, we found the essential ingredients to 
be the interaction and the restriction of y^ to Ly^ that 
is periodic boundary conditions in the ^/-direction. The 
detailed form of the interaction determines the bias cur- 
rents for which the steps occur. 
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single vortex response. The excess vortex moves in an 
even more nonlinear viscous fashion compared to / = 

0- 

Now we include an ac component to the dc driving cur- 
rent, and we find steps at = 2n{27Tu). The factor 
two can be explained using the scenario proposed in : 
first the excess vortex pushes the vortex in front of it to 
jump, and then follows suit. In total two jumps have 
been made. Small subharmonic V = n{27:v) (n = 1,3) 
steps can be distinguished in the I-V shown in the inset 
of Fig.|. 
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FIG. 9. I-V obtained by driving an excess and a missing 
vortex, separated by Ax, in the vortex lattice of an over- 
damped fully frustrated 32 x 32 array. The ac component 
has frequency = ^ and amplitude iac — 0.10. The curves 
are plotted for separations Ax = 8, 10, 12, 14. In the inset 
we show an I-V of an extra vortex moving in the ground- 
state vortex lattice of the fully frustrated 16 x 16 array. The 



frustration is / ■ 



0.10, ^, and/3c 



0. 



FIG. 8. I-V of a vort ex- ant i vortex pair with separation 
distance Ax = 7 in an overdamped 32 x 32 array. The ac 
component of the driving current has frequency ^ = ^ and 
amplitude iac = 0.10. The labels mP signify the underlying 
periodicity of the motion on these steps (i.e. the m in Eq. (5)). 
In the upper inset we show the spectral function of the voltage 
for the 6P step. In the lower inset we show the I-V obtained 
from the model for vortex- antivortex pair motion in the text. 
We used the following parameters: Ly = 32, 27rzy = 1.186, 
B = 0.0674, and Ax = 8. Here we plot, instead of V, L^uo/V. 
Steps then arise at the values ^. The higher n the smaller 
the step width. 



We find that the I-V of an excess vortex and a missing 
vortex (with respect to the checkerboard vortex lattice 
groundstate) in the f = \ case is similar to that of a vor- 
tex and an antivortex respectively in the / = case. An 
excess vortex can be annihilated by a missing vortex. We 
calculated the I-V for different values of Ax. In these I- 
V's one again finds the fractional vortex- ant ivortex-pair 
steps as in the / = case. One observes in Fig. ^ that as 
the separation Ax becomes larger, these steps disappear, 
and the steps corresponding to single vortices grow, as 
one would intuitively expect. Finally only the latter re- 
main. 



B. Defects in fully frustrated / = ^ array 



It was found in 11,0 that the I-V of a single excess 
vortex in a dc driven f = \ array resembles that of sin- 
gle vortex in the / = case up to a current idc = 0.34. 
At the latter current the entire checkerboard of vortices 
depins, producing a voltage that is much larger than the 



V. SUMMARY AND DISCUSSION 

In this paper we discussed the phase-locking behavior 
of a single vortex and a vortex- ant ivortex-pair under the 
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influence of an ac+dc drive. We obtained harmonic and 
subharmonic single- vortex voltage steps, and studied the 
microscopic vortex dynamics. We found period doubled 
vortex motion. The n = 1 step corresponds to a vor- 
tex jumping one plaquette per period, since in a periodic 
array all the plaquettes in the same column along the 
periodic direction are equivalent, and one would expect 
period one behavior. The simulation, in which period 
doubling vortex motion is observed, was performed in an 
8x8 array with /^c = 25. In this case the vortex motion 
can excite spin waves ]T3| , p!^ . These spin waves cause a 
spatially modulated environment for the vortex. Similar 
phase-locking is observed in a ring of Josephson junctions 

The spin waves may also be responsible for the compli- 
cated behavior of the hysteresis we observed in ac+dc 
driven arrays as a function of /3c. In recent experiments 
subharmonic giant Shapiro steps ^ were observed in 
over- and under damped arrays. In simulations it was 
shown that one can have such steps in a 2-D array by 
either including disorder ||] or inductive effects 0, or 
generally by any mechanism that breaks the translational 
invar iance of the array. In over damped single junctions 
one can only get subharmonics by using non-sinusoidal 
current-phase relationships [|6|. Our simulations show 
that a single vortex itself may phase-lock on subharmonic 
steps. The simulation of the pair case shows that single- 
vortex voltage steps are replaced by fractional steps, 
which can be understood in terms of vortex interaction. 
Unfortunately this insight also makes it manifest why 
generalization to the dynamics associated with Shapiro 
steps that involve many vortices (such as the axisymmet- 
ric coherent vortex states) is not feasible. Another conse- 
quence of the fact that we flnd subharmonic steps, is that 
a more realistic vortex equation of motion (than Eq. (|lo| ) 
and (|ll|)) should contain higher harmonic corrections to 
the sinusoidal potential. The presence of higher harmon- 
ics in the potential experienced by vortices was already 
noted by Lobb et al. 1^. 

We have shown that the parameters for which single- 
vortex voltage steps can be observed are within the reach 
of present day experiments. 
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period 
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TABLE L The periodicity of the voltage Vp averaged over 
one period of the driving force as a function of idc- The data is 
obtained from simulations of an 8 x 8 array, with parameters 
/3e = 25,iy = 1/50, and w = 0.10. 
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